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Abstract 

One important step in the renormalization group (RG) approach to a lattice 
sandpile model is the exact enumeration of all possible toppling processes of 
sandpile dynamics inside a cell for RG transformations. Here we propose a 
computer algorithm to carry out such exact enumeration for cells of planar 
lattices in RG approach to Bak-Tang-Wiesenfeld sandpile model [Phys. Rev. 
Lett. 59, 381 (1987)] and consider both the reduced-high RG equations pro- 
posed by Pietronero, Vespignani, and Zapperi (PVZ) [Phys. Rev. Lett. 72, 
1690 (1994)] and the real-height RG equations proposed by Ivashkevich [Phys. 
Rev. Lett. 76, 3368 (1996)]. Using this algorithm we are able to carry out 
RG transformations more quickly with large cell size, e.g. 3x3 cell for the 
square (sq) lattice in PVZ RG equations, which is the largest cell size at the 
present, and find some mistakes in a previous paper [Phys. Rev. E 51, 1711 
(1995)]. For sq and plane triangular (pt) lattices, we obtain the only attrac- 
tive fixed point for each lattice and calculate the avalanche exponent r and 
the dynamical exponent z. Our results suggest that the increase of the cell 
size in the PVZ RG transformation does not lead to more accurate results. 
The implication of such result is discussed. 
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I. INTRODUCTION 



In 1987, Bak, Tang and Wiesenfeld (BTW) proposed the concept of self-organized 
criticality (SOC) in order to understand the automatic (i.e. without a tuning parameter, 
such as temperature) appearance of abundant self-similar structures and scaling quantities 
in nature. BTW also proposed a lattice sandpile model and used Monte Carlo simulations 
to simulate this model on square and simple cubic lattices. They did observe self-similar 
structures and scaling quantities in the simulation data without tuning any parameter. Since 
1987, many natural phenomena have been related to SOC, such as earthquakes 0, forest 
fires |§, biological evolution |§, rice pile dynamics ||, turbulence ||, etc. Many lattice 
models have also been proposed to illustrate the behavior of SOC or avalanche processes 
iTjII . It has been found that the BTW's sandpile model is Abelian |9j and some quantities for 
this model could be calculated exactly |9|-|l2|. The BTW's Abelian sandpile model (ASM) 
|J has been considered to be a prototypical model for SOC. Many ideas about behavior 
of SOC models, such as universality and scaling [|n|], or new methods for studying SOC 
models, such as renormalization group theory [|14|-|T7|1, are often first tested in the BTW's 
ASM. In the present paper, we propose a computer algorithm which is useful for carrying 
out renormalization group transformations for the BTW's ASM. 

The BTW's ASM on a lattice 3? of N sites is defined as follows. Each site of 3? is 
assigned a height integer; the i-th site is assigned Z{ for 1 < i < N. In the beginning of 
the simulation, the height integer at each site is randomly chosen to be 0, 1,. . . , or z c — 1, 
where the critical height is z c — 1 and z c is the coordination number of the lattice. For each 
time interval T a , one particle falls on a randomly chosen lattice site, say the i-th site; the 
height Zi then becomes Zi + 1. If the new Zi < z c , then randomly choose a site again, say the 
k-th site, to add a particle; if Zi > z c , then the i-th site topples and its height z% becomes 
Z{ — z c . At the same time, each of the nearest neighbor (nn) sites of the i-th site receives 
one particle, i.e., — > Zw(ij) + 1, Vj where uj(i,j) is the label of the jth nn site of the 
i-th site. This relaxation procedure takes time t w and we assume t w /T a — > 0. If some of 
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the new z u ^j), V j, are equal or larger than z c again, these sites are denoted by 
Then, the toppling process continues in parallel for each j' with z^uji) — > z u uji\ — z c and 
z u{u)(i,j'),k) Zui(w(i,j'),k)+1 where j'), k) is the k-ih nn site of the j'-th nn site of the z-th 

site. The relaxation time for these parallel toppling and receiving processes between 
and u>(uj(i,f), k) takes another time t w . Usually, the open boundary conditions are used so 
that when a boundary site topples, the particle can leave the system. The dynamical process 
continues until the heights of all sites are less than z c . In general, if the last toppling site is 
u)(u(...(u)((jj(i, ii), 12), h)---), in), the total toppling process takes time n x t w . In this way, a 
series of toppling processes with toppling area s (i.e., the total number of toppled sites) and 
relaxation time nt w = t appears and forms an avalanche which has no characteristic size. 

After repeating many times the process of adding one particle on a randomly chosen 
site with subsequent relaxation when the height of the site is equal or larger then z c , we 
can obtain a distribution of toppling area P(s) and calculate the average relaxation time 
< t > for avalanches with toppling area s. It has been found that P(s) ~ s~ T with the 
avalanche exponent r and < t >~ s z l 2 with the dynamical exponent z. Manna |TBJ used 
Monte Carlo simulations to calculate r = 1.22 and z = 1.21 for the BTW model on the 
square (sq) lattice. Majumdar and Dhar conjectured that z = 5/4 which is consistent with 
their own numerical simulations |19|], and Priezzhev et al. [^Oj proposed that r = 5/4. By 



scaling argument, Tebaldi et al. |2l| suggested r = 6/5. Many investigations have been just 
focused on numerical simulations or exact results for height probabilities and exponents for 
the sq lattice and there is little attention to other kinds of lattices. It is not clear whether 
the sandpile model on two dimensional lattices have the same set of critical exponents. 

The renormalization group (RG) theory has been used successfully to calculate critical 
exponents, order parameters, etc, of ordinary phase transition models, e.g. the Ising model 



the Potts model and the percolation model p4]-p6|; it has also been used to 
understand universality and scaling in ordinary phase transition models. It is of interest to 
know whether one can also use RG theory to calculate critical exponents for lattice SOC 
models, e.g. the BTW sandpile model. In 1994, Pietronero, Vespignani, and Zapperi (PVZ) 



14| first proposed a RG theory to calculate critical exponents of the BTW sandpile model 
on the sq lattice. They used a stationary condition to provide a feedback mechanism that 
drives the system to its critical state. In fact, in their approach, the height of a site Z{ is 
reduced to only three cases corresponding to three classes: stable for Zi < z c — l, critical for 
z i — z c ~ 1) and unstable for z% > z c . The obtained renormalization group equations allow 
them to get an attractive fixed point. At this fixed point, they obtained avalanche exponent 
t =1.253 and dynamical exponent z =1.168. In 1995, PVZ |15| described in more details 



their RG approach to the BTW sandpile model. In 1996, Ivashkevich |L6[] generalized PVZ's 
RG scheme with real height parameters Zj by kinetic equations and used the technique of 
generating function to construct RG equations to obtain the exponent r =1.248 and height 
probabilities. In addition, Ivashkevich found that for the sq lattice critical particle transfer 
probabilities are close to branching probabilities for spanning tree. However, Papoyan and 
Povolotsky [|Tj]] pointed out that such connections could be wrong and used the example 
of the plane triangular (pt) lattice to illustrate their ideas. In 1996, Vespignani, Zapperi, 
and Loreto [^7] proposed a real space dynamical driven renormalization group to provide 
a theoretical basis for previous RG studies of the ASM. In 1997, Lubek and Usadel used 
extensive Monte Carlo simulations to find that r « 1.33 [^S] . 



One important step in the RG approach to a lattice sandpile model is the exact enumer- 
ation of all possible toppling processes of sandpile dynamics inside a cell for RG transforma- 
tion. PVZ and Ivashkevich divided the height configurations inside a RG cell into subsets 
with the same number of critical height by hand. After this procedure, for each height 
configuration, PVZ counted the possible toppling events by hand and Ivashkevich did this 
by the technique of generating function. Based on Refs. [|1~4]-|T6|| , in the present paper we 
formulate a computer algorithm for RG studies of the BTW sandpile model on the sq, pt, 
and honeycomb (he) lattices. The advantage of our method is that we don't need to classify 
all height configurations inside a RG cell into different subsets by hand which is different 
for different kinds of lattices and RG cell sizes. Therefore, it is easier for us to extend the 
same algorithm to different cell sizes and different kinds of lattices. Another advantage of 



our method is that we can study larger RG cells. For example, in |14|,[16],|29| and ||30|| , they 
have only carried out 2x2 cell to one site and 5-site cell to one site RG transformations for 
the sq lattice, and we can carry out 3x3 cell to one site RG transformation (RGT) for the 
sq lattice. 

In the present paper, we consider both PVZ and Ivashkevich RGTs. Using the PVZ 
reduced height RG equations, we calculate the density of critical sites p, the probability 
for one site to transfer sands to i different nearest neighbors pi, the avalanche exponent 
t and the dynamical exponent z by 2 x 2 and 3x3 cells to one site RGTs for the sq 
lattice, and 3 sites cell to one site RGT for the pt lattice. Using the Ivashkevich's real 
height RG equations, we calculate the critical height probabilities, the critical sand transfer 
probabilities, the avalanche exponent r, and the dynamical exponent z by 2 x 2 cell to one 
site RGT for the sq lattice and 3 sites cell to one site RGT for the pt lattice. We find that 
our calculated height probabilities are consistent with the exact and numerical simulations. 

The outline of this paper is as follows. In Sec. II, we introduce the reduced height RG 
equations. Here, we extend PVZ's formulation for the sq lattice to a general lattice. In Sec. 
Ill, we present a systematic computer algorithm for the exact enumeration of all possible 
toppling processes of sandpile dynamics inside a cell for RGT. By using this algorithm, we 
present the obtained fixed points of reduced height RG equations and the critical exponents 
in Sec. IV. In Sec. V, we apply real height RG equations to the BTW model on sq and 
pt lattices. For a different definition of RG cell for the he lattice, in Sec. VI we calcu- 
late the critical sand transfer probabilities and find that they are different from branching 
probabilities of spanning trees on the same (he) lattice. This result confirms Papoyan and 



Povolotsky's result JT^] for the pt lattice. Finally, we summarize our results and discuss 
problems for further studies in Sec. VII. 



5 



II. REDUCED HEIGHT RENORMALIZATION GROUP EQUATIONS 



In this section we extend PVZ's formulation of RG transformation [H] for the sq lattice 
to a general lattice 3? and present a simplified version which is useful for realization by a 
computer program. Then, we propose an algorithm for counting events of RG procedure. In 
the next section, we will show that this calculation is basically consistent with that of Ref. 
|T5| for a 2 x 2 cell of the square lattice. However, there are some mistakes in Ref. [[HJ . 

As in ordinary real space RG transformations for the spin and percolation models |2"2]-|2l)f . 
in the RG transformation for the BTW sandpile model on the lattice 3? with lattice constant 
a, we first divide 3? = S^ - 1 into g x g cells, then consider the transformation from each g x g 
cell into a cell with a smaller number of sites or a supersite. In the present paper, we consider 
only the later case. From the relationship between the cell and the supersite, the properties 
of the specified cell can be approximately represented by the supersite. After the first step 
of RG transformation, we have a lattice, call 3?' = of supersites with lattice constant 
ag. We can carry out the second step RG transformation on 3?' to obtain a new lattice 
3?" = 3?( 2 ) with lattice constant ag 2 and we can continue such RG transformation to obtain 
a series of lattices 3?^ with lattice constant ag k . Examples of RG cells for sq and pt lattices 
are shown in Fig. 1(a) and 1(b), respectively, whose linear dimensions (i.e. g values) are 2 
and a/3, respectively. 



Here we briefly review PVZ's RG approach to the BTW's ASM p^J15| . If we consider 
a sandpile RG cell with m = g 2 sites, the initial height Zj of each site i inside a RG cell 
can take integer values 0, 1, . . . , or q — 1, where q is the coordination number of lattice 3?; 
for example Figs 1(a) and 1(b) correspond to the sq and the pt lattices, with q = 4 and 6, 
respectively. After adding one particle to one lattice site, say site i, if the site i is unstable 
(i.e., Zi > z c ), the site i will topple and transfer particles to different nearest neighbors. 
In general, e« can take values 1,2,..., or q after RG transformation which is different from 
6i = q for the original sandpile evolution on As an illustration, the pictures of e« for 
the sq lattice are shown in Fig. 2. The toppling rule can be described by the transferring 
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probability p = (pi,p 2 , ■■■,Pq) where p ei is the probability that the unstable site i transfers 
particles to e$ different nearest neighbors; if is smaller than q, one should also specify which 
directions the particles go into, which will be explained in more details in the next section. 
In general, the evolution of the ASM at the initial stage $t(°\ is characterized by the initial 
height configuration set Z = {z} ={(z 1 , z 2 , z m )\ Zi — 0, 1, or q — 1, Vi} and the toppling 
rule e— (e±, e 2 , e m ) =(q, q, q). That is pi — p 2 — ... = p q -i = and p q = 1 for p at 
this stage. For example, for the ASM model on the original square lattice, the height and 
transferring probability can take value Zi = 0, 1, 2 or 3 and p = (0, 0, 0, 1), respectively. On 
the other hand, after the RG transformation, the height configuration set and toppling rule 
configuration set are Z and E, respectively, where E = {e} = {(ei, e 2 , e m )\ = 1,2,..., 
or q, Vi}. 

In this way, the transferring probability p and the initial height probability 
n=(n , ni, n q -i), where rii is the probability that the height of the supersite is i, are 
useful to characterize the properties of the coarse grained sandpile dynamics. We use p( fc ) = 
{Pi\p2 k \ ■■■■>Pq k ^) an d n* 1 ') = (n^, n q k l]) to denote the transferring probability and 

initial height probability on lattice (n( k \p^) and (n^ -1 ), p( fe_1 )) are linked by the 

k-th RG transformation. Instead of pf^ 1 = for i < q and and Y^Zlnf^ = 1 in 

the original level, the normalized condition J2i=iPi^ — 1 an d Sf=o = 1 are use< i i n the 
coarse grained level. 

In order to simplify the RG transformation, PVZ reduced the height of a site z,- L to only 
three cases corresponding to three classes: hi = which is stable for Zi < z c — 1, hi — 1 
which is critical for — z c — 1, and hi = 2 which is unstable for Zi > z c . Therefore, Z will 

— # 

be replaced by H — {h} = {(hi, h 2 , h m )\ hi = or 1, Vi}, The evolution rule is as follows: 
If a critical site % with height hi — 1 receives one particle, its height is now hi = 2 and 
site i becomes unstable. Then, this unstable site will transfer particles to neighboring sites 
with the transferring probability p. After this action, this toppled site becomes stable with 
hi = 0. In this approach, the stable site is still stable even if this site receives more than 
one particle. It means that the multi-toppling process of one site is omitted in PVZ's RG 
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approaches. Since PVZ used height configuration H and toppling rules configuration E to 
describe the RG sandpile evolution, the height probability vector n can be simplified to one 
parameter p, which is the probability of critical sites, and (n^,pW) becomes (p^ k \p^). 

An avalanche event specifies how a toppling process goes on. Figure 3 shows a typical 
avalanche event of a sandpile evolution on a three-site cell of the pt lattice. The labels of 
sites and toppling directions are shown in Fig. 3(a). Figure 3(b) shows a starting height 
configuration with site 1 being critical, site 2 being critical and site 3 being stable, i.e., 

— * 

h = (h±, h,2, h 3 ) = (l, 1, 0). This starting height configuration appears with a probability 
p 2 (l — p). When we add one particle to site 1, site 1 becomes unstable and transfers one 
particle to site 2 via directions 2 and out of cell via direction 6. Note that the probability 
of adding one particle to site 1 is 1/3 since there are 3 sites inside the RG cell and the 
probability of transferring particles to directions 2 and 6 from site 1 is P2/C2 since there 
are Cf possible toppling rules for transferring two particles to six directions. After site 2 
receives one particle from site 1 it becomes unstable and transfers one particle to site 1 via 
direction 5, one particle to site 3 via direction 6, and two particles out of cell via directions 
1 and 2; the probability of transferring particles to four nn sites from site 2 via directions 1, 
2, 5, and 6 is P4/C4. Since site 1 will not topple again and site 3 is a stable site, the toppling 
process stops. In summary, the probability of the avalanche event and the relaxation time 
(t) of Fig. 3(b) are (l/3)p 2 (l - p)(l/15)p 2 (l/15)p 4 and 2 (we set t w = 1), respectively. 

Under RG transformation, each three-site cell of -ft^ -1 ) transfers into one supersite of 
as shown in Fig. 1(b) and Fig. 3(c); a site of 3?( fe_1 ) can transfer particles to six 
directions shown in the right-hand side of Fig. 3(a) and a supersite of can transfer 
particles to six directions shown in the right-hand side of Fig. 3(c), which correspond to six 
outgoing directions of Fig. 1(d). In the avalanche event of Fig. 3(b), the particles which 
leave the cell in directions 6, 1, and 2 contribute to outgoing particles in directions 5, 1, and 
2, respectively, in the figure at the right-hand side of Fig. 3(c) (see also Fig. 1(d)). Thus 
the avalanche event of Fig. 3(b) in ^ k ^^ contributes to for the supersite of 3?^. 

In general, the probability of one avalanche event is determined by the initial height 
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configuration, the starting site of the avalanche, and the detail of the toppling process. If 
the initial height configuration is the j-th element of H and contains m! critical sites, the 
starting site is site a, the toppling process consists of a& unstable sites which transfer k 
particles to k different directions for 1 < k < q, the probability of this avalanche event is 

P E (j,a, ai ,a 2 ,...,a q ) = W,W <a> (^p 1 ) a ^^p 2 ) a \..(^ Pq ) a \ (1) 



where Wj = p m (1 



p) m m ' and W <a> = 1/m, ai,a 2 , and a q are zero or positive integer 



and Yh=i a i = m " which is the total number of toppled sites. Since we don't consider 
multiple toppling events (i.e., one site can topple at most once), m" should be equal to or 
smaller than m. 

To construct RGT, we should first collect all possible events on a RG cell of ^ k ^ 1 \ 
Consider the set of events Cj(a, ax, a 2 , . . . , a q , t) which are initiated by the j-th initial height 
configuration of H and starting site a, have relaxation time t, topple times by toppling 
rule with k directions. For the set Cj(a, ax, a 2 , . . . , a q , t), we count the number of events on 
3fj( fc_1 ) which topple i directions in point of view of site on $iS k \ This number is denoted by 
Bi(j, a, ai, a 2 , a q , t). Therefore, for 1 < i < q 



/ <q i i i \ f^qV2 ) • • • \*~fi 



qPq 



(k-l)\a q 



a ai,a,2,...,a q t 

E E Y. B itii a i a ii a 2,---,a q ,t)P E (j,a,a 1 ,a 2 ,. ■ .,a q )/W> 

a ai,a2,...,a q t 



(2) 



is the summation of probability of toppling events on K^ -1 ^ which evolve under the j-th 
initial configuration on K^ -1 -* and contribute to the on 3fjW. Note that it is possible to 
have toppling processes inside a cell but transfer nothing outside a cell. PVZ omits all of 
these events and Ivashkevich et. al. |31| simply assumed these events appear with probability 
Pq. Now the properties of the system are fully characterized by the distribution (p^ k \p^) 



at this scale. Then, the RG equation can be written as 
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(3) 



where n is the total number of configurations in H which could contribute to the right-hand 

side of Eq. (|3|), Wj, the weight of the j-th element of H, is a function of p^ -1 ), N w and Fj 

are the normalized constants of Wj and fij with N w = J2j Wj and Fj = J2i fij > respectively. 

This RG equation can be well understood by the transition rate and master equation of 

sandpile evolution P7| , |ST]] . 

Another important assumption is the inflow of particles equals the flow of particles out 

of the system, i.e., in the stationary state, < p >= [jT4]j3l]l . This implies p^' -1 ^ is a 

function of p^ for j from 1 to q and can be written as 

1 

Ei i x p\ 



P {k ~ 1] = ~ . (»-!) • (4) 



Replacing p^ -1 ) in Wj and of Eq. (|3|) by Eq. (^), p\ k ' is a pure function of pf l ' for 
j from 1 to q. From RG assumption, at k — > oo, the critical transferring probability p* = 
(p*)P2> •••>Po) an( i the critical density p* can be obtained. 



III. A SIMPLE ALGORITHM FOR COUNTING EVENTS 

In lattices for k > 1 obtained after RG transformations, one site can transfer particles 
to one, two, or q different directions with probabilities pi, p%, and p q , respectively. 
There are totally q different directions which can be labeled by di, d>2, and d q . In a 
more precise description, a site can transfer one particle or nothing to a specified direction 
di which is denoted by the variable with values = 1 or 0, respectively. Therefore, 
the toppling rules of sites are exactly determined by the vector r = (r dl , r d2 , r dq ). Since 
one unstable site sends nothing to its neighbors is omitted in PVZ RG transformation, 
(r dl ,rd 2 , r dq ) = (0, 0, 0) = is forbidden. For a given r, the site will send j = J2f=i r di 
particles to j different directions. If we divide the toppling rule configuration set S = 
{{r dl ,r d2 , ...,r dq )\r dt = or 1 , Vi and Ei=i r * ^ °} into subsets % = {(r dl ,r d2 , ...,r dq )\r di = 
or 1 , Vz and J2l=i r di — j} according to the value j, q subsets can be obtained. And it 
is straightforward to show that the number of elements of the j-th subset $Sj and 3 are CJ 
and Yfj=\ C] = 2 9 — 1, respectively. 
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Now we define Pro.(r) as the probability that an unstable site transfers particles with 
a given toppling rule r. The relation pj = J2r& 3 Pto.(y) can be obtained. On the other 
hand, we assume that the probability Pro.(v) is the same for all r G because of the 
isotropy of lattice 9ft. Therefore, we conclude that Pro.(r) = Pi/Cf for r G Q'j. For example, 
Pro.(r = (1,0,. ..,0)) = Pro.(r = (0, 1,...,0)) =...= Pro.(r = (0,0,...,1)) = pjq and 
Pro.{r = (1,1,. ..,0)) = pilC\. In this way R={(Fi, r 2 , r m )|Fj G Vi} can be used to 
represent the toppling rule configuration of an RG cell; R contains more information than 
E = {e} introduced in Sec. II. 

In Sec. II, we presented PVZ approach to construct RG transformation of Eq. (|3]) by 
using avalanche events and we used Fig. 3 as an illustrative example. In this section, we 
will use the idea of toppling rule configuration R and a computer algorithm to calculate 
all terms which contribute to the right-hand side of Eq. (|]). Consider a RG cell of m 
sites, whose initial high configuration is (hi, /i 2 , h m ) and the toppling rule configuration 
is (r!,r 2 , ...,r m ). If we add one particle to site a, the sandpile evolution starts and finally 
evolves into one state where all new hi are stable. In order to count all possible avalanche 
events, we prepare a Fortran computer program SOCRG to generate such events. In SOCRG, 
we consider three sub-configurations sets and use (m + 2) do-loops to generate such sub- 
configurations: 

1. The initial height configuration H: There are two possible states for each site i.e., 
hi = or 1. Therefore, totally, 2 m possible configurations are considered. However, 
some of these 2 m configurations will not induce sandpile evolution, for example, hi = 
for all i. In SOCRG, we use the first do-loop to generate all possible configurations in 
H. 

2. The starting point set of sandpile dynamics ^ = {a\a = 1,2,..., or m }, where a 
denotes the site on which we add one particle: There are m possible positions to add 
one particle to a RG cell with m sites. In SOCRG, we use the second do-loop to 
generate all possible site a. in ^f. If the chosen site a is not a critical site, SOCRG 
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goes to the next site. 

3. The toppling rule configuration R: There are 2 q — 1 possible toppling rules for each 
site. Therefore, there are (2 q — l) m possible toppling configurations for an m-site cell. 
In SOCRG, we use the third to the (m + 2)-th do-loops to generate such toppling 
configurations. 

For a specified configuration with h s , a s , (rf , r|, . . . , r£j generated by these (m + 2) do- 
loops, we can calculate the toppled vector O = (0±, 2 , ■ ■ ■ , O m ) where Oj=l or depending 
on whether site i topples or it doesn't during the sandpile evolution. For example, Fig. 3 
shows an event on a three-sites RG cell of the triangle lattice. In this example, (h{, h 2 , h^) 
obtained from the first do-loop, a s obtained from the second do-loop, rf obtained from the 
third do-loop, and r 2 obtained from the 4-th do-loop, are (1, 1, 0), 1, (0, 1, 0, 0, 0, 1), 
and (1, 1 ,0, 0, 1, 1), respectively. Then, O = (1,1,0) is obtained. The weight of this 
height configuration is p 2 (l — p). The probability of this event is l/3p 2 (l — p)Pro.(r{ = 
(0, 1, 0, 0, 0, l))Pro.(r 2 = (1, 1, 0, 0, 1, 1)) and the relaxation time t — 2. Finally, this event 
is classified as r = (1, 1, 0, 0, 1, 0) for the supersite of Fig. 3(c). In this case, since site 3 does 
not topple (0 3 = 0). the probability of the toppling process depends only on h s , a s , rf , and 
?2- The do-loop used to generate F3 can be skipped quickly (see below). 

In general, consider the ix-th element of the first do-loop h 11 , the ^2-th element of the 
second do-loop a 12 , and the jk-th. element of the (k + 2)-th do-loop r jk for 1 < k < m; 
the combination of such elements is denoted by the D{i 1 ,i 2 , ji, . . . , jk, ■ ■ ■ ,jm), which is a 
configuration of the set consists of H, \I>, and R. This configuration will generate a toppling 
event with probability 

m 

P P {h}\a'\{v^}) = W Kil W <ah> H(Pro.(P*))°*. (5) 

k=i 

Since Ok = means that site k is not involved in the toppling process, any choice 
of r in the (k + 2)-th loop will correspond the same event represented by the config- 
uration D(ii, iii ji, ■ ■ ■ ) jk, ■ ■ ■ , jm)- Define a set K which consists of the configurations: 

12 



D(ix, z 2 , Ji, . . . , Jk, ■ ■ ■ , J m ) with Jfc = jk for Ok = 1 or J for = 0, where J is an 
integer and 1 < J < 2 q ~ 1 . From O, we can calculate the total number of toppled 
sites: m" = 0\ + 2 + • ■ • + O m . There are (2 q — i) m_m " elements in K and every 
such element of K corresponds to the same event which is generated by the configuration 
D{i\i H,j\, ■ ■ ■ , jk, ■ ■ ■ ,jm)- In SOCRG, the do-loops corresponding to sites with Ok = are 
passed through quickly to save the computing time and the following technique is used. 

Define an m-dimensional array: V(ji, j 2 , • • • , jm), where 1 < jk < 2 9_1 for 1 < k < m. 
Immediately after a new height configuration is chosen by the first do-loop (say the i\- 
th step) and a new starting point is chosen by the second do-loop (say the 22-th step), 
all elements of the array V are set to 1. For each D(ii, i 2 , j%, . . . , jk, ■ ■ ■ , j m ) config- 
uration generated in the third to the (m + 2)-th do-loop, we first check the value of 
V(ju327 ■■■ ,jm)- If V(ji, ji, ■ ■ ■ , jm)=l, we record this event and obtain a set K. Then, 
we set V(Ji, J 2 ,...,Jm)=0 for D(ix,i 2 , Ji, J 2 , • • • , J m ) in K. If V {31,32, ■ ■ ■ , jm)=0, we skip 
this step to the next step under the third do-loop to the (m + 2)-th do-loop. In this way, the 
do-loops corresponding to sites with Ok = can be passed through quickly. Repeat above 
procedure for different combinations of height configuration and starting point and we finish 
the calculation of all do-loops. 

Basically, the form of Eq. (|j) can be transferred to the form of Eq. ([!]). Therefore, 
we can construct the RG Eqs. @ and @. In order to test this algorithm, we calculate 
RGT for a 2 x 2 cell to one site on square lattice, which is shown in Figs. 1(a) and (c) 
and has been done in details by Vespignani, Zapperi, and Pietronero (VZP) [plj. Define 
Ai{k, a t , a 2 , a 3 , a 4 ) = £ a £t B^j, a, a x , a 2 , a 3 , a 4 , ^)(^) ai (^) a2 (^) a3 (^) a4 x S{k,4>{j)) 
where 0(j) is the number of critical sites of the j-th configuration of heights. We find that 
almost all of our calculated values of Ai(k, a±, a 2 , 0,3, a 4 ) are the same as those which appear 
in the appendix of JT§, except A 3 {3, 0, 1, 1, 1) = 0.749997, A 2 (4, 0, 4, 0, 0) = 0.1772839, 
and ^(4,2,0,1,1) = 1.460250 in |T§. Our values of A 3 (3, 0, 1, 1, 1) , A 2 (4, 0, 4, 0, 0), 
and A 4 (4,2,0,l,l) are 1.750000, 0.1172839, and 1.406250, respectively. We find that 
A 2 (4, 0, 4, 0, 0) and A 4 (4, 2, 0, 1, 1) in Ref. [15| are only slightly different from our values. 
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They might be typos in ||15|| . However, there is an obvious difference between our value and 
VZP's value for A 3 (3, 0, 1, 1, 1) We believe that VZP's value is wrong since all other 240 
terms of VZP are exactly the same as ours and VZP calculated Ak(i, ai, 02, £13, 04) by hand 
and we calculate Ak(i, a\, 02, 03, 04) by a systematic algorithm. 

IV. CALCULATIONS OF CRITICAL EXPONENTS 

In the following calculations we use the RG approach in the form of Eqs. (|3]) and (f|). 
First, we consider the RG calculation from the small scale transformation for square and 
triangle lattices which are shown in Fig. 1. All bonds outgoing from a given cell into another 
are included into one renormalized bond of the supersite, as is shown in Figs. 1(c) and 1(d). 
It is one of the simplest choices where supersites of Figs 1(c) and 1(d) contain just a few 
sites of Fig 1(a) and Fig. 1(b), respectively. A block of sites in J^ -1 -* (i.e. an RG cell) 
is replaced by a site of W k \ Here g is equal to 2 and \f?> for square and triangle lattices, 



respectively. A more complex choice for the square lattice can be found in [|30| where one 
supersite contains five sites. In the present paper, for the square lattice we will extend this 
study to the case where one RG cell contains nine sites. This is the largest RG cell which 
has been considered for the RG approach to the BTW sandpile model. 

By using the algorithm of section III, we enumerate all possible toppling events for two 



kinds of cells shown in Fig. 1. However, as in |l^JT5|| we drop those events in which the 
number of toppled sites are smaller than g t =2; here it is assumed that such events contribute 
to po which has been discussed in Sec. II. In other words, these events are assumed to 
transfer no particle to nearest neighbor cells. Of course, this is an approximation. From the 
procedures stated above, we can express Wj and fij of Eq. @ in term of p( fc_1 ) and p( fc_1 ). 
The normalized factors W n and Fi can also be calculated. By using initial values of p*- ** 
and p^ to iterate Eqs. (|3|) and (^) until k — > 00, we can obtain the fixed point of p and p 
shown in Table I. Note that our results for 2x2 cell to one site RG transformation for the 
square lattice are slightly different from those of Refs. |in|,|r^J. This is due to the difference 
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in values of Aj(i, aj, 02, 03, 0,4) discussed in the Sec. III. For each lattice, we choose some 
different starting points of p^ and p(°). We find that all starting points evolve into the 
same fixed point which is shown in Fig. 4. This shows that there is only one attractive fixed 
point of RG equation for square and triangle lattices in this RG transformation. 

In WM, the critical exponent of the avalanche distribution r in P(s) ~ s~ T is calculated 
at the fixed point of the RG transformation through the fixed point parameters p* and p*. 
First, the probability that one avalanche occurs on 3t( k ~^ but doesn't occur on 8v k > can be 
written as 

# = X>*(1-PT (6) 

i=i 

In addition, the probability K can also be understood as the probability that the linear 
dimension of the avalanche, I, is larger than ag k ~ 1 ( the lattice constant of 9?( fc_1 ) ) and 
smaller than ag k (the lattice constant of 3?^). Thus P(s)ds = P(l)dl ~ s~ T ds =Z 1 ~ 2T (i/, 
where the toppling area s ~ / 2 , and we have 

K = W = 1 - o 2 ^ (7) 

also holds. Using Eqs. (H) and (0), we obtain r which is listed in the part (A) of Table II. 
We find that our r for square and triangle lattices are close to the numerical simulations 
28| , |33[| for the square lattice. 

Another independent critical exponent is the dynamical exponent z. From the scaling 
laws at the fixed point we know that the average time of a dynamical process scales with the 
linear length as < t >~ l z . Therefore, the time scale of a relaxation event on lattice ^ k ' 
and tk-i on the lattice 3?( fc-1 ) are related by the relation tk/tk-i = (agV*> / ag( k ~ x ') z = g z . On 
the other hand, the time scale can be obtained as a function of the time scale t^-i from 
the RG equations. The relation is given by % =< t^ 1 ^ > t^i, and 

^ ^ Wi W a Bj(i,a,ai } a 2 , ...,a q ,t k -i) 



< > = £ 



N w F { 



x(^r i) ) a H^?" i) ) a2 ---(4pf- i) ) a '(vi) (8) 

Vj2 
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where < t^'" 1 ) > is the average number of subprocesses on ^ k ^^ needed to have a relaxation 
process on 3?w. By inserting the fixed point parameters into the calculation of < t >= 
linifc^oo < £W >, we obtain the following result for the dynamical exponent fl4 

In < t > 

Using Eqs. (0) and <H), we obtain 2; which is listed in the part (A) of Table II. We find 
that our result z = 1.284 for the triangle lattice is not far from the theoretical prediction 
value z = 1.25 [B5|. If the universality is valid for square and triangle lattices, the obtained 
z must be the same both lattices. In our RG calculations, we find z for the square lattice 
has a larger deviation from the theoretical predicted value: z = 1.25. 

We also consider 3x3 cell for the square lattice, which is shown in Fig. 5. Here g t = 
3. The critical density of sites and transferring probabilities are shown in Table I. We find 
that there are some difference between the p* and p* in 2 x 2 cell to one site and 3x3 cell 
to one site RG transformations. The exponents of r and z shown in Table II have larger 
deviations from the simulation and theoretical prediction results than the results obtained 
from 2x2 cell. 



V. REAL HEIGHT RENORMALIZATION GROUP EQUATIONS 

In the above study, the height configuration Z is simply characterized by H. In this 
section, the real height configuration Z is used to build the RG equation. Instead of (p, 
p) in above calculations, (n, p) is used in this section, where n = (no, ni, n g _i) with 
n o + n i + ••• + n q _i = 1 and rij is the probability that the height is i. Therefore the Wi 
of Eq. (1) is replaced by Wi = irp =1 n Zr And the steady state equation of Eq. (|J) is now 
characterized by the following conditions: = h\ = ... = h q -\ = [|1^]. Therefore, the 
relationships between concentration of height n^ k ' and transferring probability p^" at the 
stationary state is: 

n£\ = (£p^-)/pW (10) 
i=i 
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where p^ = n| =1 i x p( k ' is the average number of particles sent from one site to other sites. 

By using the computer algorithm in section III, we count all events over Z ® ^ eg) R. 
Here, for height configuration (zi, z 2 , z m ), there are q possible states for each site, i.e., 
Zi = 0, 1,..., or q — 1. Totally, q m possible configurations are considered. Therefore, the 
number of events with real height Z is much larger than that with H. According to Eqs. ([3]) 
and (|TUD, the p* and n* can be obtained after repeated iterations. In Table III, we compare 
the critical height probabilities with the exact ]nj and numerical results. The numerical 
results are obtained from simulations on 1000 x 1000 sq and pt lattices and in each case 10 6 
configurations are generated to obtain the data. We also list the results of RG calculations 
reported in Refs. [|I^,|l7[]. We find that our RG fixed point is very close to previous RG 



calculations [jT6| , |lT|| . It means that our RG calculation is reliable. We also find that the RG 
results of the critical height probabilities n* have the same behavior as the numerical or 
exact results. 

For critical exponent r, Eq. (|7|) is still useful but Eq. @ should be revised to satisfy 
the definition of n. The new equation is 

K = Y,p*{l-nir (11) 

-i=i 

where p* in Eq. @ has to be replaced by n* By Eqs. (0) and (|TTD, the critical exponent r is 
obtained for sq and triangle lattice which are shown in the part (B) of Tab. II. We find that 
the value of r is very close to previous calculation for sq and triangle lattice. It shows again 
that our computer algorithm and calculation are equivalent to Ivashkevich's algorithm. In 
Refs. fll6yi7f , they didn't calculate the dynamical exponent z. Here, by Eqs. (H) and (|j), we 



calculate this quality which is shown in part (B) of Tab. II. We find that the obtained z on 
triangle lattice is far from the value of exact and numerical results for the square lattice. 

VI. CALCULATIONS FOR THE HONEYCOMB LATTICE 



There is no overlap of sites between RG cells in the Figs.l and 5. In other words, 
from these two figures, if one site belongs to one specified cell, this site does not belong to 
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other cells. Due to this property of cells, Eqs. (|q)-(|Tl|) can be used to calculate critical 
exponents. Consider the RG transformation cell for the honeycomb (he) lattice shown in 
Fig. 6. This kind of RG cell allows one site to belong to two different RG cells. If we use 
the cell transformation in Fig. 6(a) and 6(b), we can still obtain the fixed point of p and pi. 
However, it seems that it is inappropriate to get the critical exponents r and z for the he 
lattic by Eqs. (|)-(0). 



The algorithm in Sec. Ill is still useful for the calculation of transferring probabilities 
for the honeycomb lattice. First, for the reduced height RG equation, we count all possible 
toppling events and drop those events with g t = 3. Then, using Eqs. @ and (|J) to obtain 
p* and p* which are shown in Table I. For real height RG equation, the obtained n* and 
p* are shown in Table III and IV, respectively. We find that the height probabilities are 
quite consistent with the numerical simulations on 1000 x 1000 honeycomb lattice, where 
10 6 configurations are generated to obtain numerical data. 

The branching probability of a spanning tree || pk is the probability for any site of a 
random spanning tree to have a coordination number k. The Green functions of the Laplace 
equation for the square, triangle, and honeycomb lattices are well known [M\. For an infinite 



lattice, G depends only on the difference of site coordinates rather than their values. Here 
we show the Green function for the honeycomb lattice as an example: 

[Tl T2) Jo Jo 21- UAcos 2 a4:Cosacos 8) 27r27r' 1 ' 



where 



2 

f(x 1 ,x 2 ,yi,y 2 ,a,P) = cos(a(xi - x 2 )P(yi - y 2 ))-cos(a)cos(a(x 1 - x 2 ) 

(3(yi - y 2 ))^cos(a(x 1 - x 2 )P{yi - 3/2)/?) - 1- (13) 



Following the same procedure as Ref. |P^j , we obtain the branching probability of a spanning 
tree pk for the honeycomb lattice. The values of pk for sq and pt lattices |T7|,|52"| are also 
listed in Table IV. In Table IV, we compare the RG transferring probabilities p* with the 
branching probabilities of spanning tree pi. We find that they have the same behavior on 



the square lattice. However, they are not on the triangle and honeycomb lattices. Thus, we 
can conclude that the hypothesis about the coincidence of Pi and pi proposed in [[[(J is not 
valid. 



VII. SUMMARY AND DISCUSSION 

In this paper we use a computer algorithm to calculate the effective toppling events for 
two kinds of RG equations. We find that the values of critical exponents r and z for 2x2 
cell to one site RG transformation for the sq lattice are closer to conjectured exact values 
and simulation values than those for 3 x 3 cell to one site RG transformation. There are two 
possible reasons. The first reason is that we don't consider the multiple toppling events in 
the RG calculations. The errors arising from the multiple toppling events are larger for 3x3 
cell. Therefore, when the cell size increases, we should consider the multiple toppling events 
in order to get accurate values of r and z. This can not be done easily in exact enumeration 
approach to RG transformations, but can be done by Monte Carlo RG calculations. Another 
possible reason is that the critical exponent r is not well defined for the BTW sandpile model 
and therefore RGT can not be used to obtain the critical exponent. 

For the PVZ approach, the critical density of sites, p*, is equivalent to height probabil- 
ities, n* in the Ivashkevich approach. If we compare the value of p* with the numerical or 
exact n*, we find that p* is larger than n*. This is related to the rule that in PVZ approach 
the stable site will not topple even the stable site receive more than one particles. Therefore, 
the obtained p* must be larger to compensate the loss of stable sites which have the potential 
of toppling. Again, when the cell size is larger, the compensation effect is also larger. For 
example, p* obtained from 3x3 cell of the square lattice is larger than that obtained from 
the 2x2. There is no such kind of problem in Ivashkevich's approach. We find the n* is 
consistent with the numerical simulations and exact results. 

In summary, it is worthwhile to consider the multiple toppling events in the real height 
RG treatment in order to answer the question discussed above. However, it is hard to carry 
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out exact enumerations for larger cell sizes. In the next step, we plan to use Monte Carlo 
simulations to construct RG transformations with large cells and include multiple toppling 
events. 
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FIGURES 

FIG. 1. (a), and (b). Transformation from a cell to a site for the RG transformation. Here (a) 
and (b) are for square lattice with 2x2 cells and for the plane triangular lattice with y/3 x y/3 
cells, respectively. Blocks of sites of the coarse grained lattices Jft^ -1 ) become supersites sites on 
the renormalized lattice 5?^. (c). and (d). We show that the directions outgoing from the blocks, 
which are encircled by rectangle corner roundness, are coupled to the directions of the lattice at 
the next scale. 

FIG. 2. The picture of toppling rules for a site of the square lattice. The ei denotes one 
unstable site i will transfer particles to e« different nearest neighbors. For each e*, there are CJ. 
possible toppling rules. 

FIG. 3. A toppling process on a three-site cell of the plane triangular lattice. Open dots 
represent stable sites, black dots represent critical sites, and encircled black dots represent unstable 
site, (a) The label of sites inside a cell and the label of toppling directions of an unstable site, (b) 
A series of toppling process from left to right. For site 1 and 2, their toppling rule ?i and F2 are 
(0, 1, 0, 0, 0, 1) and (1, 1 ,0, 0, 1, 1), respectively, (c) After the RG transformation, the three-site 
cell of (a) the toppling processes of (b) are represented a supersite and the toppling rule r = (1, 1, 
0, 0, 1, 0) on the supersite. Note that the toppling directions in (c) can be obtained from (b) by 
rotating 30° counterclockwise. 

FIG. 4. Iteration results for different initial values of p*- ^ and p(°\ The solid line with symbol 
o, □, A, and v represent the 2x2 cell transformation on square lattice with initial values p^ ^ = 
(0, 0, 0, 1), (0, 0, 0, 1), (|, \, \, |), and (|, \, \, \) and p(°) = 0.1, 0.8, 0.1, and 0.8, respectively. 
The dashed line with symbol o, □, A, and V represent the 3 sites cell transformation on triangle 
lattice with initial values p(°) = (0, 0, 0, 0, 0, 1), (0, 0, 0, 0, 0, 1), (±, \, ±, ±, ±, ±), and and 
(g> h h h h g) = °- 8 ' and °- 8 ' respectively, (a). A plot of the against the 
A;th iteration number k. (b). A plot of the |p( fc )| 2 = (p^) 2 + ip^) 2 + ■■■ + ipq^) 2 against the fcth 
iteration number k. 
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FIG. 5. Transformation from a cell with 9 sites on square lattice, (a). This shows the trans- 
formation from a cell to a site. (b). We show that the directions outgoing from the blocks, which 
are encircled by rectangle corner roundness, are coupled to the directions of the lattice at the next 
scale. 

FIG. 6. Transformation from a cell with 6 sites on honeycomb lattice (a). This shows a 
transformation from a cell to a site. (b). We show that the directions outgoing from the blocks, 
which are encircled by rectangle corner roundness, are coupled to the directions of the lattice at 
the next scale. 
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TABLES 



TABLE I. Critical density and transferring probabilities at the RG fixed point. 





P* 


Pi 


* 

P2 


* 

P3 


Pa 


* 

P5 


Pi 


sq [2 x 2]( a ) 


0.515 


0.327 


0.437 


0.205 


0.031 






sq [2 x 2]W 


0.468 


0.240 


0.442 


0.261 


0.057 






sq [3 x 3]( fl ) 


0.663 


0.572 


0.353 


0.070 


0.005 






pt [3 sites] ^ 


0.214 


0.000002 


0.0005 


0.040 


0.314 


0.582 


0.073 


he [6 sites] 


0.763 


0.702 


0.285 


0.013 









(a) : this work based on PVZ approach 

(b) : PVZ 1994 M 
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TABLE II. Avalanche exponent r and dynamical exponent z for square and plane triangular 
lattices. 





t (so) 


z (so) 


r (pt) 


z (pt) 


sq [2 x 21 W (A) 


1.243 


1.147 






sq [2 x 2]( 6 ) (A) 


1.253 


1.168 






sq [5 sites] ( c ) (A) 


1.235 


1.236 






sq [3 x 3]( fl ) (A) 


1.122 


1.082 






pt [3 sites] W (A) 






1.363 


1.284 


sq [2 x 2]( a ) (5) 


1.248 


1.150 






sq [2 x 2]M) (B) 


1.248 











pt [3 sites] ( a ) (B) 






1.367 


1.433 


pt [3 sites] ( £ ) (B) 






1.367 




Simulation 


1.33W 


1.254® 






Prediction 


1.25(f) and 1.2W 


1.25 W 







(A) 
(B) 

(a) : 

(b) : 

(c) : 

(d) : 

(e) : 

(f) : 



: reduced parameter RG calculations 
: real parameter RG calculations 
this work 
PVZ 1994 



Moreno, Zapperi, and Pacheco 1999 |30| 
Ivashkevich 1996 
Papoyan and Povolotsky 1997 [171 
Lubek and Usadel 1997 |2S| 
Priezzhev, Ktitarev, and Ivashkevich 1996 [20] 
Tebaldi, Menech, and Stella 1999 fH 



(h) 

(i): Majumdar and Dhar 1992 |L9 



27 



TABLE III. Comparisons of critical height probabilities n* for square (sq), plane triangular 
(pt), and honeycomb (he) lattices obtained by RG transformation, numerical simulations, and 
exact calculation. In the numerical simulations with statistics of 10 6 configurations on 1000 x 1000 
lattices are generated to obtain the data. 





* 

n 


* 

n\ 


n*. 


* 


* 


* 


RG (sq)( a ) 


0.021 


0.134 


0.349 


0.496 






RG (sq)^ 


0.021 


0.134 


0.349 


0.496 






Simulation (sq)( a ) 


0.074 


0.174 


0.306 


0.446 








Exact (sq)( c ) 


0.074 


0.174 


0.306 


0.446 






RG (pt)( a ) 


0.036 


0.135 


0.198 


0.210 


0.211 


0.211 


RG (pt)W 


0.036 


0.135 


0.198 


0.210 


0.211 


0.211 


Simulation (pt)( a ) 


0.058 


0.094 


0.139 


0.188 


0.240 


0.281 


RG (hc)( fl ) 


0.014 


0.308 


0.678 








Simulation (hc)( a ) 


0.083 


0.293 


0.624 









(a): this work based on real height RG approach 



(b) : Ivashkevich 1996 p] 

(c) : Priezzhev 1994 |HJ 



(d): Papoyan and Povolotsky 1997 [17 
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TABLE IV. 


Fixed point of transferring 


probability p* 


and branching 


probabilities 


(BP) of 


spanning trees p 
















Pi (Pi) 


P*2 fe) 


P*3 (P3) 


Pi (P4) 


Pi (P5) 


Pt (pe) 


RG (sq)( a ) 


0.295 


0.435 


0.229 


0.041 





— 


RG (sq)( & ) 


0.295 


0.435 


0.229 


0.041 





— 


BP (sq)( c ) 


0.295 


0.447 


0.222 


0.036 





— 


RG (pt)( fl ) 


0.0000179 


0.00226 


0.0558 


0.296 


0.471 


0.174 


RG (pt)W 


n nnnm 7Q 


n nn99fi 

U.UUiiU 


0.0558 


0.296 


0.471 


0.174 


BP (pt)W 


0.322 


0.417 


0.207 


0.049 


0.006 


0.0002 


RG (hc)( fl ) 


0.546 


0.432 


0.022 








BP (hc)( a ) 


0.25 


0.5 


0.25 









(a) : this work based on real height RG approach 

(b) : Ivashkevich 1996 @ 



(c) : Manna, Dhar, and Majumdar 1992 [|| 

(d) : Papoyan and Povolotsky 1997 JL7 
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